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1. Introduction 



The program of topological quantum computation is to realize fault tolerant 
quantum computation using topological phases of quantum systems [FKLW] . The 
central open question is whether or not there exist such physical systems which 
are capable of performing universal quantum computation. In [F], a family of 
Hamiltonians Ho i is proposed as candidates for the Chern-Simons phases, which 
C*~) ' are known to support universal quantum computation for each level I > 3, 1 =/= 4 

[FLW1,FLW2]. Freedman conjectures that the perturbed ground states of H 0i i 
are given by the Drinfeld double of the S'0(3)-Witten-Chern-Simons topological 
quantum held theories (TQFTs). The approach in [F] is an algebraic study of the 
i-C ■ effect of perturbation based on a rigidity result of the picture TQFTs [FNWW] . The 

picture TQFTs are the Drinfeld double of the SO(3)-Witten-Chern-Simons TQFTs. 
^ r"| The idea is to treat local relations in picture TQFTs as perturbations. In this 

paper we will investigate numerically the perturbed ground states of Hq^i for some 
■ computationally tractable cases. While numerical study of Hamiltonians is carried 

out routinely in physics literature, we face a dilemma here. The major motivation of 
our investigation is quantum computing, but one motivation of quantum computing 
as Feynman pointed out is to study quantum systems numerically. Therefore, we are 
in a self-referential situation. This is also manifested in the fact that our numerical 
computation quickly reaches the limit of present computing power. 

Freedman's Hamiltonians H$ i define quantum loop gas models on any celluated 
compact surface. We study the simplest nontrivial cases: celluations of the torus. 
Our numerical data support Freedman's conjecture, but the conjectured space of 
ground states does not come out in full. Study of substantially larger systems is 
necessary, but computationally intractable now. One new phenomenon we discov- 
ered is some lonely states in the ground states of i?o,(- Those lonely states give rise 
to unwanted ground state vectors of Hqj which persist in the perturbed ground 
states. There are several possible explanations of those phenomena: the small size 
of the system, the choice of our perturbation, or the Euclidean geometry of the 
torus. We also observe clearly the expected energy gap between ground states and 
the first excited states. 

The Hamiltonians Hq i are 7 th order interaction Hamiltonians. It is an open 
question to find 2 th order or 3 th order interaction Hamiltonians with approximately 
the same ground states. Note that [F] also contains a family of 4 th order interaction 
Hamiltonians similar to Hq,i- 

i 
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2. Freedman's Hamiltonian 

At the heart of this study lie the Hamiltonians Hq,i, for levels I > 1, which grants 
the structure of a TQFT to the hypothetical physical systems. Each Hamiltonian 
H ,i is described as a sum of local projections, which "implement" the concept of 
combinatorial generalized isotopy {g(d) -isotopy) on a celluated surface. 

2.1. Combinatorial Isotopy. Fix a closed oriented surface S, and let A be a 

triangulation of E with n vertices. The dual graph A* to A therefore defines a 
celluation of £ by n polygons, with the edges of A* being boundaries of adjacent 
cells. 

A spin configuration is an assignment s : A* ^ {+, — }ofa positive or negative 
spin to each cell in A*. We denote a "spin flip" by ~: on entire spin configurations, 
s denotes a global interchange of + and — ; for a cell c of A*, s c denotes the spin 
configuration which agrees with s away from c, and flips the spin at c. 

A spin configuration can be thought of as a 2-coloring of S, partitioning £ into 
+ and — regions. We will be interested in studying the boundary d s of these + 
and — regions with respect to spin configurations s. We will refer to d s as domain 
walls of a spin configuration. Since A* is dual to a triangulation, the domain walls 
are all 1-manifolds in X. Our goal will be to capture global properties of domain 
walls d s in terms of local manipulations at each cell. 

Thus, for a fixed cell c, we will need to consider the boundary dc of c, which 
consists of edges in A*. Color each edge on dc according to the neighboring cell; 
set d + c to be the edges which bound a +-ccll, and 9_c to be the edges which 
bound a — cell. We then say that the pair (s, c) is type-g if both d+c and d-c form 
a connected topological arc. Note that neither d+c nor d-c can be empty. We 
define the pair (s, c) to be type-h if c and all neighboring cells have the same spin 
(d s(c) c = dc). 

The motivation for these definitions is as follows: if (s, c) is type-g, then the +/— 
boundary d s meets dc along one of the topological arcs d + c or d-c (the one which 
is opposite the spin of c). By changing the sign of c, we replace the part of d s which 
meets dc with the other of the topological arcs. But this transformation can be 
viewed as a fixed endpoint isotopy of one arc to the other, with isotopy occurring 
over the 2-cell c (see the figure below). So, if (s,c) is type-g, then d s and d S c are 
isotopic. 




For type-h pairs (s,c), the cell c and all of its neighbors have the same spin. 
So, if the spin is flipped on cell c, then a closed loop is added to the boundary d s . 
Type-h cells are therefore responsible for the "generalized" aspect of g(d)-isotopy, 
namely, the insertion or removal of disc-bounding loops. 



ON FREEDMAN'S LATTICE MODELS FOR TOPOLOGICAL PHASES 



3 



Using this terminology, we can consider a combinatorial form of g(d)-isotopy, 
relative to the celluation A*. Given a spin configuration, we can apply: (i) type-g 
moves, which consist of flipping the spin of a type-g cell; and (ii) type-h moves, 
which consist of flipping the spin of a type-h cell. We say that two spin configura- 
tions are combinatorially g(d)-isotopic if one can be reached from the other via a 
sequence of such moves. Then we have the following proposition: 

Proposition 1. Let s and t be combinatorially g(d)-isotopic spin configurations 
with respect to the celluation A*. Then d s is g(d)-isotopic to dt- □ 

The comments above show that this result is immediate. The converse is not 
quite true: it is possible to have non-combinatorially g(d)-isotopic spin configura- 
tions s and t such that d s and dt are g(d)-isotopic. The reason for this is that the 
tiling could be too coarse to allow room for g(d)-isotopy to take place; for example, 
on a hexagonal tiling of the torus, a spin configuration which assigned different col- 
ors to adjacent vertical rows would have no type-g or type-h cells (but its boundary 
would be g(d)-isotopic to that of its dual). 

2.2. Definition of Hoj. We are now ready to consider the definition of the Hamil- 
tonian Hq i. For a triangulation A with n vertices, we associate the 2™-dimensional 
Hilbcrt space H = ®" =1 C 2 ; let ci,...,c n denote the dual 2-cells in A*. Then 
we can express a basis for TL by {|s}}, where s runs over spin configurations and 
\s) = |s(ci)) <8> |s(c 2 )} <8> • • • <S> |s(c„)}. Set the parameter d = 2 cos j^, and put: 

H , l= ]T \s-s c )(s-s c \ + ]T l s -^" C )( s -^" C | 

( S ,C) (s,c) 

typo g typo h 

We interpret these terms as follows: H i is "indifferent" to g(d)-isotopy, hence the 
equal weighting on the type-g terms. The factor of ^ expresses that "loops are 
worth d" , as flipping the type-h configuration creates a closed bounding loop in the 
spin configuration. 

The primary observation, to be proved in the next section, is that the ground 
state Go,; of H ,i consists of equivalence classes of combinatorially g(d)-isotopic 
spin configurations. 

3. Theoretical Analysis 

First we consider analytic results relevant to the conjecture. The ground state 
of H ,i is determined explicitly. 

3.1. The Ground State Go,;- For a spin configuration s, let n(s) denote the 
number of trivial (disc-bounding) closed loops in s. Thus if s' is the configuration 
s with all trivial closed loops removed, we have d n ^s' ~ s. 

Theorem 2. Fix a surface E, a triangulation A of 12, and let ~a be the equivlance 
relation of combinatorial g-isotopy on spin configurations with respect to A. Then 
the ground state Go,; of the corresponding Hamiltonian H ,i has as a basis: 

K)= £ d^\ s ) 

where d = 2 cos an d the s\, . . . , s k are representatives of the k equivalence classes 
determined by ~a- 
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Proof. First, observe that the subspace V, = span{|s) : s ~a «i} is invariant under 
Hoj; this is immediate from the definition of Hq.i, as each projector acts invariantly 
on one of the Vi and trivially on Vj for i ^ j. So we can write TL = ® 4 Vi, with H ,i 
acting invariantly on each summand; therefore, the ground state can be likewise 
decomposed, Go,; = ® i ker(i7 ,z|v;). So the theorem reduces to the claim that 
ker(£/o,zk) = span^). 

First, we show that \vi) lies in ker(i?o,z|v;). For any s ~a Si, we have: 

(s\H ,i\vi)= d n(s )-d n ^+ rf " (s) - \ dn[SC) + 

cells c s.t. cells c s.t. 

(s, c) type g (s, c) type h 

cells c s.t. cells c s.t. 

(s c ,c)typog (S",c) type h 

For any (s, c) that is type-g, we have n(s) = n(s c ), since the two configurations 
are (strictly) isotopic. Similarly, if (s, c) is of type-h, then n(s c ) = n(s) + 1, so that 

d n(s) _ l d n(s") = Q. when (g C;C ) ig u typc h „ _ thcn ^ d «( S ) _ l rf «(g c ) = Sq &U 

terms in the summation cancel; since this holds for all s ~a Si, this shows that 
\vi) E kcr(i7 ,z|v t )- 

The converse claim follows similarly; for, suppose some nonzero vector \v) E Vi 
is in the kernel of i7 ,i|vi- Then (s\v) ^ for some s ~a Si, so that H j\s}(s\v} 
will have nonzero components along each \s') which differs from s by one type- 
g or type-h move. So in order for \v) € ker(i?o.;|v , i ), it must also have nonzero 
components along these \s). Continuing this argument, we see that (s\v) ^ for 
all s ~a Si- Analyzing the above computation, we see that if s and s' differ by a 
type-g move, then (s\v) = {s'\v}; and if they differ by a type-h move (say, s' is s 
with a loop removed), then d(s\v) = (s'\v). These requirements in turn force \v) to 
be a multiple of \vi), so that ker(i/ ,z|vj = span |t>j). □ 



4. Numerical Study 

Celluations of a surface S dual to triangulations have a nice topological property. 
The domain walls in a spin configuration arc all 1-manifolds. The homology class 
represented by those domain walls is the zero class as the domain walls are bounding 
1-manifolds. The linear combinations of domain walls coming from different spin 
configurations in a fixed celluation forms a finite dimensional vector space. Those 
vector spaces are combinatorial approximations of a picture TQFT. Local relations 
in picture TQFTs are the topological realization of perturbations. The rigidity of 
the picture TQFTs says the only non-trivial local relations for each I is generated 
by the Jones-Wenzl projector. Therefore, for each level I perturbations of the 
Hamiltonian i?o.z will result in only one possibility: the picture TQFTs at level 
L The purpose of the present study then is to perform numerical analysis on 
computationally tractable cases, in order to test this conjecture. We will focus on 
perturbations of the form H e = H + eV. As a preliminary case, we will choose 
V to the sum of a x on each cell. Our study focuses on / = 3, so the loop value 
is d = 2 1 ■ This corresponds to the double of the 50(3)-Witten-Chern-Simons 
theory at level=3. In [FLW1], we have shown that the level I = 3 theory for both 
SU{2) and SO(3) supports universal quantum computation. 
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hex 18a 



Figure 1 . The hexagonal tilings of the torus used in the numerical study. 

4.1. On Tilings. The numerical study focused on tilings of the torus. Nine hexag- 
onal tilings of the torus were considered, representing the "obvious" hexagonal 
tilings of the torus which are computationally tractable. These tilings are depicted 
in Figure^ on the torus visualized in the standard way as a rectangle with opposite 
sides identified. 

The tiling hex7 is the dual of a minimal triangulation of the torus, and thus is a 
minimal tiling. The other tilings are, loosely speaking, formed by p adjacent verti- 
cal columns of q cells; when p is odd, an extra "twist" is needed to align the vertices 
properly. Thus, tilings with (p, g) = (3, 3), (4, 3), (3, 4), (5, 3), (3, 5), (4, 4), (3, 6), (6, 3) 
are represented. 

4.2. Ground State Vectors of H i. First, a combinatorial computation was 
performed to verify the result of Theorem|21 regarding the form of the ground state 
vectors of H Q i. Spin configurations were grouped into g-isotopy classes, and vectors 
corresponding to each class were created according to the formula of Theorem [21 
Each such vector was found to be in the ground state of Hqj; further, the number 
of such (orthogonal) vectors was equal to the dimension of the ground state, as 
calculated via an eigenvalue computation, so that the ground state was numerically 
verified to be exactly the span of such vectors, in all cases considered. 

The interesting result of these calculations concerns what we will call lonely 
configurations: spin configurations which are not g-isotopic to any other spin con- 
figuration. (In other words, configurations in which there are no type-g or type-h 
cells.) Some of the tilings considered admit such lonely configurations, and others 
do not; in considering the numerical analysis of the perturbed Hamiltonians below, 
it will be necessary to take these lonely configurations into account. 

We considered nine different hexagonal tilings of the torus, ranging from 7 cells 
(the minimal hexagonal tiling of the torus, corresponding to the minimal triangu- 
lation of the torus) to 18 cells (the maximum computationally tractable case). The 
following table presents the results of this preliminary analysis of the tilings: 

4.3. Perturbation Ground States G e j. For the numerical analysis of the per- 
turbed ground states, we calculated the lowest energy eigenvalues of H e ^, for e 
ranging from to 1 . Plots of the eigenvalues as a function of e are given below, for 
each of the nine tilings in the table above. 

As the collected data simply provides the lowest eigenvalues, the plotted lines 
simply indicate the trajectories of the various eigenvalues. The lines themselves do 
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Table 1. Ground states of H j, with respect to given tilings 



Tiling 


n 


dim G , 3 


Lonely configurations 


Non Lonely 


hcx7 


7 


5 





5 


hex9 


9 


5 





5 


hcxl2a 


12 


8 


2 


6 


hcxl2b 


12 


17 


12 


5 


hexl5a 


15 


7 





7 


hcxl5b 


15 


8 





8 


hexl6 


16 


24 


18 


6 


hexl8a 


18 


16 


8 


8 


hcx!8b 


18 


21 


14 


7 



not follow a particular eigenvalue, but rather, they connect the lowest eigenvalue, 
the second- lowest eigenvalue, etc.; in other words, the plots below should be treated 
more like scatterplots. The connecting lines are included to help elucidate the 
trajectories of the eigenvalues and make it easier to determine where eigenvalues 
converge and diverge. 

Also, the plots do not indicate the number of eigenvalues represented by a par- 
ticular line. For this, the numerical data had to be examined by hand. Thus all 
lines represent one eigenvalue, except where explicitly labelled otherwise. Not all 
multiplicities in higher energy eigenstates are labelled. 

Finally, the numerical algorithm for calculating eigenvalues did not always con- 
verge, and thus provided possibly erroneous results for particular values of e. In 
most cases these were isolated, and therefore can be safely ignored; however, there 
are a few larger regions in which the algorithm failed to converge. These also are 
indicated in the plots below, by regions surrounded by dashed lines. 

First, we consider the tiling hex7: 

Tiling hex7 




.1 I 1 1 1 1 1 

0.2 0.4 0.6 0.8 1 



The initial 5-dimensional ground state splits into four seperate states, with the 
second lowest of these being doubly degenerate, as indicated in the figure. 
For hex9: 
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Tiling hex9 




As in hex7, the initial 5-dimensional ground state splits, although this time into 
five seperate states. Again only the lowest of these remains in the ground state as 
e grows. (The indicated region of failed convergence is for 0.01 < e < 0.02.) 

The spikes appearing in the plot below for hexl2a arc due to failed convergence 
of the algorithm; however, the algorithm only failed to converge for some of the 
higher energy eigenvalues. Therefore it is not indicated as a potential error region 
(as we are interested mainly in the lowest energy eigenstates) . 



Tiling hex12a 




Though it is difficult to see, six of the eight original ground state eigenvalues 
split in the initial region (e < 0.01), with the two lowest staying together. Then 
the two lower ones fracture apart around e w 0.2. Recall also that this tiling has 
two lonely configurations, although it is not possible to tell from these data how 
the lonely configurations might correspond to the eigenvalue trajectories. 

The tiling hexl2b has 12 lonely configurations; thus the plot below indicates 
more clearly how these relate to the other ground states. 

From this plot it is almost certain that the 12 eigenvalues that initially are 
negative correspond to the 12 lonely configurations of the tiling. 
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Tiling hex12b 




Neither of the 15-tilings admit any lonely configurations. The plot for hexl5a 
follows: All seven of the ground state eigenvalues split, with only one of them 



Tiling hex15a 




2 



0.2 0.4 0.6 0.8 1 



remaining in the ground state for e > 0.3. The area of failed convergence is 0.09 < 
e < 0.3. 

For hexl5b: 

Again, all eight of the ground state eigenvalues split, and only one remains in 
the ground state beyond e w 0.25. Though it is difficult to sec on the graph, the 
line labelled with multiplicity two is not originally in the ground state when e = 0. 

The plots for the 16- and 18-tilings are restricted to < e < 0.5, as they are 
computationally expensive, and the trajectories of the eigenvalues seem clear. For 
hexl6, the general trend regarding the lonely configurations remains consistent; 
interestingly, the other pattern seems to hold, although in a slighty more dramatic 
fashion. The eigenvalue corresponding to the lowest non-lonely configuration first 
rises sharply, then hits a peak around e w 0.07, and then begins to decline, re- 
maining below the other non-ground state trajectories until e w 0.3. Although 
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Tiling hex15b 




Tiling hex16 

0.5 r : 1 1 1 1 1 




it appears to nearly "merge" with the non-ground state trajectories at e « 0.07, 
manual inspection of the data indicates that in fact it remains clearly distinct. 

The plot for hexl8a follows: Again, eight of the ground state eigenvalues sep- 
arate distinctly from the rest, almost certainly corresponding to the eight lonely 
configurations. And again, only one of the original ground state eigenvalues remains 
consistently low, with only the lonely configurations below it. 

Finally, the tiling hexl8b demonstrates the same pattern. The fourteen dimen- 
sions of the ground state corresponding to the lonely configurations initially split 
off and remain consistently lower than the rest, and one of the remaining ground 
state configurations seperates and remains in a lower state than everything but the 
lonely configurations. 

5. Discussion 

• The expected dimension of the perturbed ground states for I = 3 is 4. 
Ignoring the lonely states represented by lonely configurations, the un- 
perturbed ground states form a space of dimensions 5,6,7, and 8. There 
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Tiling hex18a 




0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 



Tiling hex18b 




0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 



was consistently one branch of the original ground state which split off from 
the rest and remained in the ground state. This is an indication of what 
we're looking for, but not coming out in full because of the small size of 
the tilings. Only two of the tilings (hexl8a and hexl5b) even admitted 
four parallel essential circles, and even then only in one direction (e.g., just 
horizontally, not vertically). 
• The problem of lonely configurations might be only a technical issue. If 
they exist, then they will certainly be the lowest energy eigenstates under 
perturbation. It seems likely the lonely configurations are "isolated" in the 
sense that, if the system is prepared in some suitable initial state, then it 
will never "jump" to one of these lonely configurations, but rather slowly 
fall into one of the ground state configurations corresponding to a large 
g-isotopy class. 

Given a triangulation of a surface, we can always subdivide the surface 
to get a new triangulation without any lonely configurations. The problem 
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with this solution is that we have to introduce much more tilings into the 
problem, which is definitely computationally intractable. 
• Recall that in the Chern-Simons theory, the Lagrange has a quadratic term 
AAdA and a cubic term |^4 3 - If the cubic term is treated as a perturbation, 
then it is a 3rd order term. So maybe a better perturbation should be at 
least of order 3. 

6. Gap in thermodynamic limit 

Inspecting the data, we observe that the lowest eigenvalue and the next one has 
a clear gap. The question is whether this gap persists if we go to finer and finer 
tilings. Further study will be carried out. 

Acknowledgement: Research of Z.W. is supported by NSF grant CISE/EIA- 
0130388 and Army Research Office. 
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Appendix A. Implementation Details 

The open source numerical package ARPACK [7] was the main tool used for 
finding the ground state and low energy eigenstates of the Hamiltonians. ARPACK 
is optimized to find certain eigenvalues (e.g., those of the lowest magnitude) of large, 
sparse, symmetric, real- valued matrices, and is thus well-suited to our problem. A 
set of templates interface the Fortran ARPACK code with C++ was used; all of 
the custom code for our problem was done in CH — h 

For a fixed Hamiltonian, let n denote the number of tiles in the associated tiling; 
thus, the dimension of the associated Ho will be 2™. 

The sparseness of the Hamiltonian can be analyzed by determining the number 
of type-g and type-h spin configuration pairs associated with a given tiling. It 
is relatively easy to determine that, for a hexagonal tiling with n tiles, there are 
62 • n ■ 2™~ 7 type-g and type-h spin-configuration pairs, and therefore twice as 
many nonzero entries in the matrix (total of 2 2n entries). While this is relatively 
sparse, it is still inefficient to compute and store the matrix explicitly (using the 
"Compressed Sparse Column" format of ARPACK), as the number of required 
entries grows exponentially. Computing the matrix-vector product "on the fly" is 
much more efficient. This algorithm is given in pseudo-code below: 
//helper method 

Bar(c,t); // returns the configuration c with tile t flipped 
H0_product (v [] ,w [] ) { //computes H_0|v>=|w> 
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w=(0,0 0) ; 

foreach tile t { 

foreach type-h cfg (c,t) { 

w[c]+=v[c]-(l/d)v[Bar(c,t)] ; 

w[Bar(c,t)]+=(-l/d)v[c]+(l/d"2)v[Bar(c,t)] ; 

} 

foreach type-g cfg (c,t) { 
w[c]+=v[c]-v[Bar(c,t)] ; 
w[Bar(c,t)]+=v[Bar(c,t)]-v[c] ; 

} 

} 

} 

He_product (v [] ,w [] ) { //computes H_\epsilon I v>= I w> 
H0_product(v,w) ; //first do the H_0 part 
foreach cfg c { // then add the perturbation 
foreach tile t { 

w[Bar(c,t)]+=v[c] *epsilon; 

} 

} 

} 

The algorithm requires 0(n2 n ) floating-point operations (flops) per matrix-vector 
product; compared with the standard product using the CSC-stored format, which 
would required 0(2 2n+1 ) flops. The configurations that are type-h and type-g, 
with respect to each tile, are calculated beforehand and stored. Therefore, the 
storage requirements are 62n2™ -7 integers, which is on the same order as the storage 
requirements for storing the matrix in CSC format. 

An attempt was made to parallelize the code, as there is a parallel version of 
ARPACK available. However, the code only allows for parallelization of the matrix- 
vector product; and the lag in communicating results across nodes would be signif- 
icantly more costly than the actual computation of the results. Thus the only way 
in which parallelization might be useful is if the actual eigenvalue algorithm was 
parallelized. This approach was not pursued. 



